Well use an utility script with few lines of code to parse the output of GRASS commands and make new functions that use the parsed output.

The script use ipython specific syntax like !system_command which allows to run any command available in the user $PATH. The code is saved in a file with .ipy extension and is imported using the ipython magic function %run:

%run file.ipy

the source code available in grassutil.ipy in the current directory and can be loaded in the notebook using the magic function %load:

%load grassutil.ipy

It is also possible to save a new file or to overwrite an existing one, from the content of a code cell using the magic function:

%%file filename.

Running the grassutil.py code, the following functions will be immediatly available in the current notebook:

getLayerList
list2dict
vlayerInfo
rlayerInfo
region2dict
makeImage

In the IPython notebook is possinle to execute any command available in the $USER $PATH environment, such all the grass command after exporting the grass environment.

Use of the g.gisenv


In [38]:
!g.gisenv


GISDBASE=/home/main/notebooks/data/grass7data
LOCATION_NAME=nc_basic_spm_grass7
MAPSET=user1
LOCATION=/home/main/notebooks/data/grass7data/nc_basic_spm_grass7

Use of the g.mapset

  • Set the GRASS WORKSPACE to :
    • LOCATION : spearfish
    • MAPSET : user1

In [39]:
!g.mapset location=nc_basic_spm_grass7 mapset=user1


WARNING: <user1> is already the current mapset

  • print projection info with g.proj

In [40]:
!g.proj -p


-PROJ_INFO-------------------------------------------------
name       : Lambert Conformal Conic
proj       : lcc
datum      : nad83
a          : 6378137.0
es         : 0.006694380022900787
lat_1      : 36.16666666666666
lat_2      : 34.33333333333334
lat_0      : 33.75
lon_0      : -79
x_0        : 609601.22
y_0        : 0
no_defs    : defined
-PROJ_EPSG-------------------------------------------------
epsg       : 3358
-PROJ_UNITS------------------------------------------------
unit       : Meter
units      : Meters
meters     : 1
  • list vector and raster layers with g.list

In [41]:
!g.list rast


basin_50K
basins
elevation@PERMANENT
elevation@user1
elevation_shade
geology
lakes
landuse
soils
  • use the getLayerList function to store the g.list output in a python list

In [42]:
rasterlist = getLayerList(type='rast')
vectorlist = getLayerList(type='vect')

In [43]:
rasterlist


Out[43]:
['basin_50K',
 'basins',
 'elevation@PERMANENT',
 'elevation@user1',
 'elevation_shade',
 'geology',
 'lakes',
 'landuse',
 'soils']

In [44]:
vectorlist


Out[44]:
['boundary_region',
 'boundary_state',
 'census',
 'elev_points',
 'firestations',
 'geology',
 'geonames',
 'hospitals',
 'points_of_interest',
 'railroads',
 'roadsmajor',
 'schools',
 'streams',
 'streets',
 'zipcodes']
  • print info for a raster layer with r.info

In [45]:
!r.info elevation@PERMANENT


 +----------------------------------------------------------------------------+
 | Map:      elevation@PERMANENT            Date: Tue Nov  7 01:09:51 2006    |
 | Mapset:   PERMANENT                      Login of Creator: helena          |
 | Location: nc_basic_spm_grass7                                              |
 | DataBase: /home/main/notebooks/data/grass7data                             |
 | Title:    South-West Wake county: Elevation NED 10m ( elev_ned10m )        |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 255             |
 |   Data Type:    FCELL                                                      |
 |   Rows:         1350                                                       |
 |   Columns:      1500                                                       |
 |   Total Cells:  2025000                                                    |
 |        Projection: Lambert Conformal Conic                                 |
 |            N:     228500    S:     215000   Res:    10                     |
 |            E:     645000    W:     630000   Res:    10                     |
 |   Range of data:    min = 55.57879  max = 156.3299                         |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.proj                                                     |
 |                                                                            |
 |   Comments:                                                                |
 |    r.proj input="ned03arcsec" location="northcarolina_latlong" mapset="\   |
 |    helena" output="elev_ned10m" method="cubic" resolution=10               |
 |                                                                            |
 +----------------------------------------------------------------------------+

  • use the 'r/v'layerInfo function to store the r.info / v.info output in a python dictionary

In [46]:
rasterlayerinfo = rlayerInfo(map='elevation')
vectorlayerinfo = vlayerInfo(map='geology')


mapset parameter not specified, using mapset PERMANENT as default
['north=228500', 'south=215000', 'east=645000', 'west=630000', 'nsres=10', 'ewres=10', 'rows=1350', 'cols=1500', 'cells=2025000', 'datatype=FCELL', 'ncats=255']

In [47]:
rasterlayerinfo.keys()


Out[47]:
['rows',
 'north',
 'datatype',
 'west',
 'cells',
 'cols',
 'range',
 'ewres',
 'ncats',
 'east',
 'nsres',
 'south',
 'history']

In [48]:
rlayerInfo('elevation')


mapset parameter not specified, using mapset PERMANENT as default
['north=228500', 'south=215000', 'east=645000', 'west=630000', 'nsres=10', 'ewres=10', 'rows=1350', 'cols=1500', 'cells=2025000', 'datatype=FCELL', 'ncats=255']
Out[48]:
{'cells': '2025000',
 'cols': '1500',
 'datatype': 'FCELL',
 'east': '645000',
 'ewres': '10',
 'history': 'Data Source:\n   \n   \nData Description:\n   generated by r.proj\nComments:\n   r.proj input="ned03arcsec" location="northcarolina_latlong" mapset="\\\n   helena" output="elev_ned10m" method="cubic" resolution=10',
 'ncats': '255',
 'north': '228500',
 'nsres': '10',
 'range': {'max': '156.3299', 'min': '55.57879'},
 'rows': '1350',
 'south': '215000',
 'west': '630000'}
  • use of the makeImage function to display raster and/or vector maps (a wrapper around display commands in GRASS using the CAIRO DRIVER)

In [ ]:
!g.mapset location=nc_basic_spm_grass7 mapset=user1
inputlayer={
    'raster': ['elevation'], 
    'vector':['points_of_interest']
}

makeImage(basemap='elevation', inputlayer=inputlayer, maptype='overlay', 
          vsize=10, maptitle='points_of_interest', gridsize=4000, outputimagename='test.png')

In [ ]:
from IPython.core.display import Image

Example on how to repoject raster and vector data between 2 different GRASS LOCATION:

  • set the GRASS environment to the previously generated lonlat LOCATION (see GRASS_init.ipynb)
  • reproject GRASS layers :
    • raster : elevation
    • vector : bugsites

In [ ]:
!g.proj -c epsg=4326 location=lonlat

In [ ]:
!g.mapset -c location=lonlat mapset=PERMANENT

In [ ]:
region = !r.proj input=elevation location=nc_basic_spm_grass7 -g

In [ ]:
region

In [ ]:
newregion = dict([(i.split('=')[0],i.split('=')[1]) for i in region[-1].split()])

In [ ]:
!g.region -p n={newregion['n']}  s={newregion['s']}  e={newregion['e']}  w={newregion['w']} res=0.0001

In [ ]:
!r.proj input=elevation location=nc_basic_spm_grass7 output=elevation method=bicubic --o --q

In [ ]:
!v.proj input=points_of_interest location=nc_basic_spm_grass7 output=points_of_interest --o --q

In [ ]:
#!g.region -p n={newregion['n']}  s={newregion['s']}  e={newregion['e']}  w={newregion['w']} res=0.0001
inputlayer={
    'raster': ['elevation'], 
    'vector':['points_of_interest']
}

makeImage(basemap='elevation', inputlayer=inputlayer, maptype='overlay', 
          vsize=10, maptitle='points_of_interest', outputimagename='test.png')

In [ ]: